Synchronization of two interacting populations of oscillators 
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We analyze synchronization between two interacting populations of different phase oscillators. 
For the important case of asymmetric coupling functions, we find a much richer dynamical behavior 
compared to that of symmetrically coupled populations of identical oscillators (IJ . It includes three 
types of bistabilities, higher order entrainment and the existence of states with unusual stability 
properties. All possible routes to synchronization of the populations are presented and some stability 
boundaries are obtained analytically. The impact of these findings for neuroscience is discussed. 



In its original sense synchronization describes the mu- 
tual adjustment of frequencies between two interacting 
oscillators jJi] • This route to synchronization differs from 
that taking place in large communities of oscillators [2| . 
Motivated by this fact, Kuramoto extended the notion 
of synchronization to a statistical theory of oscillator en- 
sembles y . The natural diversity among the cornponents 
was considered through either unimodal |j, |5|, |a| or bi- 
modal 4, 7] frequency distributions. 

In a pioneering work, Okuda and Kuramoto [l| ana- 
lyzed two symmetrically coupled populations of identi- 
cal phase oscillators under the influence of noise. How- 
ever, the important problem of synchronization between 
ensembles of oscillators remains almost unexplored, al- 
though communities of natural oscillators are usually 
composed of interacting subpopulations |3- For in- 
stance, it has been shown experimentally that synchro- 
nization arises between different neighboring visual cor- 
tex columns and also between different cortical areas, 
where synchronization processes are of crucial impor- 
tance (S|- Thus, it is a challenge to understand quantita- 
tively the routes to synchronization among macroscopic 
ensembles of oscillators. 

In this letter we study two populations of phase oscil- 
lators interacting asymmetrically, as it is likely to occur 
due to a number of reasons (asymmetric couplings, dif- 
ferent population sizes, time delays...). In addition, we 
consider the oscillators within each population to be non- 
identical. The system under study is then 
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where « = !,.. .,7V^1. Here, 9\ ' ' describes the phase 
of the ith oscillator in population 1 or 2, respectively. The 
populations are coupled internally with coupling strength 
Kp, whereas the interpopulation coupling is determined 
by K. The asymmetry is introduced into the model 
through the phase shift < a < tt/2 |J, \^ . This permits 
the oscillators to synchronize to a frequency that deviates 
from the simple average of their natural frequencies. This 



behavior is common to many living systems such as mam- 
malian intestine and heart cells |3| . Moreover, such asym- 
metry appears in the phase reduction of nonisochronous 
oscillators |3jbi| and Josephson junction arrays jlllll2l|. 
and it is used for modeling time delays l2| and informa- 
tion concerning synaptic connections in a neural network 
[l3j. The natural frequencies w- ' are considered to be 
distributed according to a density g^^''^'{uj) of width 7, 
symmetric about the mean uj^^''^' and unimodal. 

The phase coherence within each population is de- 
scribed by the complex order parameters R^^'^'e 
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which permit to write the system ^ in 



terms of the mean field quantities R^^'^'' and ■0'^^'^) 
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If the populations are uncoupled, i.e. K = 0, each 
of them reduces to the well known Kuramoto model jj] . 
For a given Kp this model exhibits a phase transition at 
a critical value of the frequency dispersal 7^. For 7 > 7c 
the oscillators rotate with their natural frequencies and 
i^(i'2) ^ 0{^yl/N), but for 7 < 7c mutual entrainment 
occurs among a small fraction of oscillators giving rise 
to a finite value of the order parameter R^^'^\ Thus, a 
cluster of locked oscillators emerges through a Hopf bi- 
furcation of frequency Q^^'^'' that, in general (a ^ 0), 
depends on the overall shape of ^(^'^^(w) |9|. The drift- 
ing oscillators arrange in a stationary distribution that 
does not contribute to the order parameters 4]. Finally, 
for identical oscillators 7 = 0, each population fully syn- 
chronizes in-phase, R^^''^'> = 1, for arbitrary small Kp. 

When K > the two locked clusters begin to interact. 
If this interaction is similar to the frequency adjustment 
between two coupled oscillators, one expects mutual lock- 
ing between these two clusters to occur via a saddle-node 
bifurcation at some K — Kc 2']. Especially, for 7 = 0, 



synchronization should arise at [Alu 
Kc = Auj/{2cosa). 
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In the following we investigate the dynamics of (O 
in the full (K, 7)-parameter plane. In the thermo- 
dynamic limit a density function can be defined so 



that p^^''^\9,t,u!)du;d9 describes the number of oscilla- 
tors with natural frequencies in [w, lo + du] and phase 
in [9, 9 + d9] at time t. For fixed uj the distribution 
p^^''^\9,t^uj) of the phases 9 is normalized to unity. The 
evolution of p^^''^^9,t,u)) obeys the continuity equation 
ap(i'2)/9i = -c»(p(i'2)^(i,2))/5g,^ fQj. ^Y^^^Yi the incoher- 
ent state po = (271')^^ is always a trivial solution |5(. The 
function p^^'-^\9,t,uj) is real and 27r-periodic in 9 and 
therefore admits the Fourier expansion p'-^-^'>{9,t,Lu) — 

EOO (1.2)/, N iW 1 (1,2) *(1,2) rrni 

i=-ooPi '{t.^y" , where p'_i ' = Pi' • ' . Thus 

the order parameter can be written in terms of the 
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Fourier components as i?^^'^^ 
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use (/(^'•^^(cl')) to denote frequency average weighted 
with (7'^^'^)([.i;), respectively). Now, after inserting jSJl into 
the continuity equation, we obtain an infinite system of 
integro-differential equations for the Fourier modes 
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The stability of po can be analyzed by studying the 
evolution of a perturbed state ^'-^'^■'(6', t, w) close to po 
(note that the p\ ' are then small quantities). Lin- 
earization of 101) reveals that the only potentially un- 
stable modes are Z = ±1 and hence I = 1 has solution 
5(1-2) (w)e^* ^ 0{\pi\^). This leads to 



FIG. 1: {K,'y) phase diagram of system Q assuming 
Lorentzian frequency distributions, Aa; = 0.5, Kp = 1, and 
for different values of a (in rad.): (a) a = 0, (b) a = 0.1, 
(c) a = 7r/4, (d) a = 1.15. Numerical stability boundaries 
(A^ = 1000) are indicated as solid lines. Dotted lines represent 
analytical stability boundaries 7c± obtained from 0. Note 
that 7c+ fully overlaps with numerical results. Region I: syn- 
chronization. Region II: coexistence. Region III: incoherence. 
Insets: bistability in the dashed regions around A (see text). 
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Considering the distribution of frequencies to be of 
Lorentzian type, f/(i'2)(w) = (7/7r)[72 + (u — w(i'2))2j-i^ 
the self-consistent problem (O can be solved analytically. 
The stability of the incoherent state po is then described 
by two pairs of complex conjugate eigenvalues, namely 
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where Z±{t) = e ('*'±)*, and c.c. denotes the complex 
conjugate of the preceding term. The modulus of the 
number zq = (Acj — VAw^ — e^'^"' K^)e~'^°' / K is a weight 
for the fraction of frequencies ri+ and ri_ in populations 
1 and 2, respectively. 

The symmetric case, a = iFiaUVa)] .— Our eigen- 
modes (fTj) coincide with those of [l| (replacing 7 by the 
noise intensity). From ® the state III can become un- 
stable in two different ways. When K > Aw the tran- 
sition III-I takes place through a single Hopf bifurca- 
tion and both populations synchronize to the same fre- 
quency CI = uj. The presence of a single macroscopic 
oscillation is denoted as region I. When K < Aoj the 
instability is through a degenerated Hopf bifurcation. 
Both (A±,Aj_) cross simultaneously the imaginary axis 

at 7c± = 7c = Kp/2 (line CA') and two macroscopic os- 
)/2 for mode I = 1, and the complex dilations with frequencies n± = =f1/2%/Acj2 - if2 ^ ^ 

emerge (note that a saddle-node bifurcation should take 
place at A', i.e. Kc = Auj). The region of coexistence of 
two different macroscopic fields is labeled as II. The inset 
of Fig^a) shows how the saddle- node line bad crosses the 
two Hopf lines 7c at a, and joins the Hopf curve 7c+ at d. 
Thus the Hopf bifurcation is supercritical all along 7c+ , 
except in aA'd, where it is subcritical and hence a region 
of bistability between states III/I and II/I is observed, 
c.c.-l- The coupling-modified frequencies of the individual os- 

cillators a)j = limt^oo 9\ ' /t, provide a useful mea- 
sure of synchronization: when iiT = (7 < 7c) the 
frequency-locked oscillators in each population form a 
(7) single plateau that is the only contribution to the or- 
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with uj = (a;(i) 

conjugate for I = —1. Imposing Re(A±) ~ defines 
explicitly the two critical curves 7c± (K) (see FigQJ. Each 
curve represents a Hopf bifurcation with frequency given 
by fl± = -Im(A±). The curve max(7c+,7c-) = 7c+ 
separates the region where the incoherent solution HI is 
stable from the unstable regions I and II. 

The eigenmodes (^p^^''^'{9,t,uj)') near criticality are 
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FIG. 2: Coupling-modified frequencies uji of populations 1 
(black) and 2 (grey) as a function of the oscillator's in- 
dex i: oscillator i has natural frequency lo^ = uj^^''^' + 
7tan[(7r/2)(2i-7V-l)/(7V + l)] (Lorentzian) {Auj = 0.5, 
id = 0, Kp ^ 1 and A^ = 1000). First row; a = 0, 7 = 0.4 
and (a) K = 0, (h) K ^ 0.4, (c) K = 0.41. Second and third 
rows: a = 7r/4, K = 0.53 and (d) 7 = 0.5, (e) 7 = 0.47, (f) 
7 = 0.187, (g) 7 = 0.15, (h) 7 = 0.12, (i) 7 = 0.118. 



der parameters (note that |zo| = in |(7J), and hence 
Tp^^') = n_t and V'^2) ^ fi^t) [Fig[2i;a)]. By increasing X, 
some of the oscillators in populations 1 and 2 begin to 
lock in a second plateau at Q-i- and f2_, respectively, ac- 
cording to (TJ [Fig|2Ib)]. Hence, i?*^^-^' begin to oscillate 
with beating frequency Afi = Q^ — ri+. Interestingly 
new clusters synchronized to higher frequencies appear 

.,(1,2) 



among those drifting oscillators with w^ ' close to 
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f7_ + nAfl, where 



1,±2,±3. 



(8) 



These plateaus fl„, which are similar to Shapiro steps 
0, 13, are not explained by Q, but they can be un- 
derstood from the fact that the drifting oscillators are 
simultaneously forced by the two order parameters iQ. 
The plateaus grow in size and in number as AQ — + 
and hence they make a nonzero contribution to the or- 
der parameters that becomes important as the system 
approaches the saddle-node bifurcating line Ba from the 
region II. At this line the synchronized state I is reached, 
Ar2 = 0, and the steps abruptly disappear [Fig|2Ic)]. 

The asymmetric case, a > [Figs^(h,c,d)]— . As a is 
increased from zero, the bifurcating lines ^c+ and 7c- 
split due to the breaking of symmetry. Interestingly, the 
eigenmodes (0 do not reflect the asymmetry through 
the amplitudes |zo|i but only through the different ex- 
ponential growths Z±{t). Figs|2Id-i) show the uJi for 
a — 7r/4 [Fig^c)], keeping K constant and decreas- 
ing 7 continuously from region III. We find that inco- 
herence [FigI2Id)] always goes unstable through a single 
Hopf bifurcation 7c+ (at Q.+ in Figl^Je)) and nucleation 
first takes place mainly within population 2. The sec- 




FIG. 3: Order parameters R-^' (black) and _R"' (grey) and 
phase difference At/j as a function of time (A'' = 1000, Auj — 
0.5, (I; = 0, Kp = 1). At i = the phases were equally 
spaced in (0, 27r]. (a) a = 7r/4, 7 = 0.12, K = 0.53 (region II) 
[FiglHh)]; (b) and (c) represent regions IF (a = 1.2, K = 0.8) 
and F (a — 1.2, K — 1.05), respectively (7 = 0, see Fig^J. 



ond Hopf bifurcation CA (at 0_ in Fig|21[f)) follows 7c- 
when the system is close enough to the incoherent state 
III. As 7 is decreased further, the system approaches the 
saddle-node bifurcation, Bd, and an increasing number 
of oscillators in population 1 becomes entrained to the 
frequencies (jHJ [Figs[2Ig,h)]. In consequence the order 
parameter R^^' oscillates with a very large amplitude 
at frequency Afi whereas i?*^^^ remains almost constant 
[FigEl^a)]. The phase difference between the order pa- 
rameters A?/) = ip^^^ — ip^"^^ reveals the presence of such 
clusters: Ai/; [Figl^Ia)] is bounded despite the fact that 
the populations are not locked in frequency [FigEfh)]. 

The bistability regions [Figsn^b,c) inset] are located 
around the intersection a of the Hopf line CA with the 
saddle-node line Bd. Within the region enclosed by Aba 
the states I and II coexist, as in the a = case. In 
contrast, the region enclosed by Aad is surrounded only 
by the state I and a new bistability between a small/large 
amplitude of the synchronized oscillation is observed. 

With increase in a, the synchronization regions I and II 
become gradually smaller because as a — *■ 7r/2 synchro- 
nization is increasingly inhibited due to frustration |2|. 
At the same time, \zo\ decreases indicating a lower de- 
gree of synchronization between the populations. This is 
in qualitative agreement with the approach of the saddle- 
node line Bd to the 7 = axis [Figs. n^c,d)]. At the crit- 
ical value a — a* the line Bb collides with the 7 = axis 
and disappears (see Inset ^d)). Therefore, for a > a* 
synchronization between the macroscopic oscillations oc- 
curs generally when the oscillation of frequency fJ- dies 
in the Hopf bifurcation CA". 

The limit 7 = [Fig[^.— The transition point B fol- 




FIG. 4; Phase diagram {K, a) for 7 = 0, Aui = 0.5. Bound- 
aries Kc and K^ are obtained analytically from Eqs.J^J and 
Q, respectively, whereas the symbols D and A correspond to 
numerical results. All other boundaries are determined nu- 
merically. Regions I (synchronization) and II (drift) are char- 
acterized by W^'^' — 1. Within regions V {Aip bounded) and 
ir {Atp not bounded) i?'^^ < 1 whereas R'-^^ = 1 (see text). 
Dashed regions present bistability between states I and IF 
(horizontal dashes) and between I and F (vertical dashes). 



lows Eq.Q as far as a < a* [Fig^. Since the oscillators 
within each population are identical, they synchronize in- 
phase, i?(^'^^ = 1, and the population's dynamics reduce 
to that of a system of two nonidentical oscillators. How- 
ever, for a > a* the synchronization transition occurs 
via a Hopf bifurcation (line A"), and the behavior in each 
population is of higher complexity. As soon as a reaches 
the critical value a* (point P), the curve Kc splits into 
two bifurcating lines, K[ and K^^ , that enclose the new 
regions IF and I' where the order parameters are not syn- 
chronized [Fig|2Ib)] and synchronized [Fig|3fc)], respec- 
tively. Within those regions the oscillators in population 

1 are not in-phase synchronized, whereas the population 

2 shows perfect in-phase entrainment [Figs|3Ib,c)]. 
Finally we outline a linear stability analysis of the syn- 
chronized state I (7 = 0) which confirms the loss of sta- 
bility of the in-phase state of population 1. From Q, 
the phase difference between the populations is Aip — 
arcsin[Aa;/(2ii'cosa)]. Then linearization of (Q results 
in a matrix with A^ — 1 eigenvalues //+ and iV — 1 eigen- 
values /i_ characterizing the stability of the in-phase syn- 
chronized state of populations 1 and 2, respectively 



^J■± 



-Kp cos a 



Kcos{±Aip + a) < 0, 



(9) 



and two eigenvalues fig — 0, ^c = — 2A'cosacos A"*/; 
(note that fic = leads to Q). Since tt/2 > Aip > 0, 
the condition is only violated for the population 1. 
Thus ^+ = determines the boundary K^ (and hence 
the point P) in very good agreement to numerics [Fig^. 
Notice that with K = we recover the in-phase sta- 
bility condition for a single population, Kp cos a > 0. 
For |a| > 7r/2 this state becomes unstable and reaches 
a neutrally-stable incoherent state. This issue has been 



the subject of a great deal of research in connection to 
the dynamics of devices consisting of Josephson junctions 
\ll\ . In the present case, however, even for |a| < n/2 the 
in-phase state in one population can be destabilized (pop- 
ulation 1) or overstabilized (population 2) due to the in- 
teraction with the other population. The global stability 
properties of the states I' and IF in population 1 are in- 
teresting directions of further study: We stress that i?^^^ 
in Figs.EIb,c) strongly depends on initial conditions and 
on perturbations, in contrast to Atp and i?/^'. 

The mean field model |^ shows rich behavior despite 
of its simplicity, especially for a ^ 0. Beyond its impor- 
tance for the theory of synchronization, oscillatory net- 
works consisting of interacting subpopulations are com- 
mon in neuroscience, and in general in many natural sys- 
tems 13- For example, synchronization seems to be a cen- 
tral mechanism for neuronal information processing and 
for communication between different brain areas 8_, 13i] . 
This plays a crucial role in the pattern recognition and 
motor control tasks Iqj. In addition, the recordings of 
neuronal activity are usually taken in different brain re- 
gions, which constitute a network of interacting subpop- 
ulations of neurons '8i|. Thus, our study represents an 
important step into understanding macroscopic synchro- 
nization in complex network architectures. 
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